function y = chirpm(fs, t, f, method)
% CHIRPM  Swept-frequency periodic waveform generator iMproved.
% adq@XJTU, last modified Aug 13 2021
% Generates swept-frequency sine/square/sawtooth/user-defined waveform. Phase
% is wrapped onto [0,2π] before any evaluation. t and f assigns  frequencies,
% fs determines Waveform base frequency, and method specifies waveform shape.
% 
% t is positive and incremental, and f is positive and of the same length.
% 
% Method should be a cell as follows:
% {'sin'/'cos'} for sin/cos;
% {'square'/'saw', duty} for square/sawtooth with specified duty cycle, where
% duty falls into [0,1];
% {'handle', @fun} for user-defined (non-numeric) 2π-periodic function;
% {'array', funValArr} for numeric function evenlly sampled over its one peroid.
% This is done with interp1 linear interpolation.
%
% The function is realized with phase-accumulating method, where ω = dφ/dt. And
% φ = 2π*int f(t)dt;
% φ(0) = 0;
% φ'(ti) = 2π*f(ti);
% Since ω(t) = 2pi*f(t) has been specified over points ti, an interpolation can
% be made, known as Ω(t); and a numeric integration yields Φ(t) that meets all
% constraints listed above.
% Different interpolation methods yield different degreeds of precision. Linear
% interp with enough points specified is desirable where function is complex.

	if fs <= 0, error('Expect positive fsamp'), end
	% if ~isscalar(method{2}) || method{2} > 1 || method{2} < 0
	% 	error('Expect scalar duty cycle on [0,1]');
	% end

	tSpan = t(end);

	tVect = 0:1/fs:tSpan;

	Phi = cumsum(2*pi.*interp1(t, f, tVect, 'linear'))./fs;
	% Phi(0) = 0 may not be met here
	Phi = wrap(Phi);

	if isequal(method{1}, 'sin')
		y = sin(Phi);
	elseif isequal(method{1}, 'cos')
		y = cos(Phi);
	elseif isequal(method{1}, 'square')
		y = Square(Phi, method{2});
	elseif isequal(method{1}, 'saw')
		y = Saw(Phi, method{2});
	elseif isequal(method{1}, 'handle')
		% isfunhdl
		y = method{2}(Phi);
	elseif isequal(method{1}, 'array')
		if ~isvector(method{2}), error('Expect vector'), end
		y = interp1(method{2}, Phi./(2*pi).*length(method{2}));
	else
		error(sprintf('Unknown parameter "%s"', method{1}));
	end

end

%---------------------------------------------------------------------------

function phi = wrap(phi)
% WRAP Remapping radiants over [0,2π)
	phi = mod(phi, 2*pi);
end

function phi = varwrap(phi)
% VARWRAP Remapping radiants over [π,π)
	phi = mod(phi, 2*pi);
	idx = phi > pi;
	phi(idx) = phi(idx) - 2*pi;
end

%---------------------------------------------------------------------------

function y = Square(phi, duty)
% Square square wave defined on [0,2π) with peaks of ±1
% 'duty' specifies the position of falling edge, [0,1]
	PI2 = 2*pi;
	y = -ones(size(phi));
	y(phi > duty*PI2) = 1;
end

function y = Saw(phi, duty)
% Saw sawtooth wave defined on [0,2π) with peaks of ±1
% 'duty' specifies the position of falling tipping point, [0,1]
% duty=0.5 means a standard triangle wave
	PI2 = 2*pi;
	idx = phi <= duty*PI2;
	y = -zeros(size(phi));
	y(idx) = phi(idx)./(duty*pi) - 1;
	y(~idx) = (PI2-phi(~idx))./(pi*(1-duty)) - 1;
end